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Abstract 

Numerical simulations of strongly correlated electron systems suffer from 
the notorious fermion sign problem which has prevented progress in under- 
standing if systems like the Hubbard model display high-temperature super- 
conductivity. Here we show how the fermion sign problem can be solved 
completely with meron-cluster methods in a large class of models of strongly 
correlated electron systems, some of which are in the extended Hubbard model 
family and show s-wave superconductivity. In these models we also find that 
on-site repulsion can even coexist with a weak chemical potential without in- 
troducing sign problems. We argue that since these models can be simulated 
efficiently using cluster algorithms they are ideal for studying many of the 
interesting phenomena in strongly correlated electron systems. 



* Current address Universiry of Utah, Salt Lake City, USA 



1 Introduction 



Strongly correlated electrons are among the most interesting condensed matter sys- 
tems. In particular, the doping of anti-ferromagnetic cuprate layers leads to high- 
temperature superconductivity. Understanding the dynamics of strongly correlated 
electron systems is very challenging because perturbative analytic techniques fail for 
systems of many strongly coupled degrees of freedom. On the other hand, numerical 
simulation methods are applicable in these cases. Unfortunately, numerical simu- 
lations of models involving fermions are far from trivial. A standard technique is 
to linearize four-fermion interactions by the introduction of bosonic fields Then 
the fermions are integrated out, leaving behind a fermion determinant which acts 
as a non-local effective action for the bosonic fields. In some cases, for example, in 
undoped cuprate layers at half-filling, the fermion determinant is positive and can 
be interpreted as a probability for the bosonic field configurations. Then standard 
importance sampling techniques can be applied to the bosonic theory. Still, due 
to the non-locality of the effective action, such simulations are very time consum- 
ing. In other cases of physical interest — in particular for doped cuprates away 
from half-filling — the fermion determinant may become negative and can hence 
not be interpreted as a probability for the bosonic configurations. When the sign 
of the fermion determinant is included in measured observables, the fluctuations 
in the sign give rise to dramatic cancellations. As a consequence, the relative sta- 
tistical errors of observables grow exponentially with the volume and the inverse 
temperature of the system. This makes it impossible in practice to simulate large 
systems in the low-temperature limit. In particular, the fermion sign problem pre- 
vents numerical simulations of the Hubbard model away from half-filling and is the 
major stumbling block against understanding high-temperature superconductivity 
with numerical methods. 

Dealing with the fermion sign problem after the fermions have been integrated 
out is incredibly complicated since the resulting bosonic effective action is non-local. 
Before integrating them out, fermions are usually described by anti-commuting 
Grassmann variables which are practically impossible to simulate directly. Here we 
propose a different strategy. Instead of working with Grassmann coherent states, we 
formulate the fermionic path integral in the Fock state basis of the Hilbert space. 
Then a fermion configuration is described in terms of bosonic occupation numbers 
which define fermion world-lines. In this case, the bosonic occupation numbers in- 
teract locally. However, besides the positive bosonic Boltzmann weight, there is a 
fermion permutation sign that results from the Pauli principle. Two fermions that 
interchange their positions during the Euclidean time-evolution give rise to a minus- 
sign. In general, the fermion world-lines which are periodic in Euclidean time define 
a permutation of fermion positions. The sign of this permutation is the fermion 
sign. The resulting sign problem is often more severe than the one associated with 
the fermion determinant. For example, in the world-line formulation a sign problem 
arises for the Hubbard model even at half-filling, while the fermion determinant is 
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positive in that case. 

Fortunately, in the Fock state basis the sign problem can be completely elimi- 
nated, at least in a large (but still restricted) class of models. This was first demon- 
strated for a system of spinless lattice fermions in @, [| by re-writing the partition 
function in terms of the statistical mechanics of closed loops with non-negative 
Boltzmann weights. Loops with a certain topology, referred to as merons, do not 
contribute to the partition function. A local algorithm was constructed in the loop 
space which avoided meron-clusters and thus solved the sign problem. It is well 
known that spinless fermions can be converted to relativistic staggered fermions by 
introducing carefully chosen phase factors with every fermion hop ||. Using this 
strategy, the solution to the fermion sign problem with spinless fermions has been 
exploited to extensively study the critical behavior of a second order chiral phase 
transition in the universality class of the Ising model || || 

Here we extend these successful techniques to fermions with spin and, in partic- 
ular, to systems in the Hubbard model family. Using a specific example, we show 
how one can build models so that the sign problem is completely solved. For tech- 
nical reasons related to the details of the cancellation of signs, the method does 
not apply to the standard Hubbard model Hamiltonian. Additional terms must be 
included. Still, the modifications do not affect the symmetries of the problem and 
we believe that the physics of some of the models discussed here is qualitatively 
similar to the standard Hubbard model. It should be pointed out that there is no 
reason to concentrate on the standard Hubbard model except for the simple form of 
its Hamiltonian. The fact that for similar Hamiltonians our method can completely 
eliminate the sign problem and yield a very efficient fermion algorithm — which is 
not the case for the standard Hubbard model — is reason enough to replace the 
standard Hubbard model by a modified Hamiltonian. After all, we want to focus on 
understanding the general physical phenomena underlying high-temperature super- 
conductivity — not necessarily on the details of one particular model Hamiltonian. 
Interestingly, in the modified model the sign problem can be eliminated even in the 
presence of an on-site repulsion at a finite chemical potential as long as a certain 
inequality is satisfied. 

1.1 The Fermion Sign Problem 

Let us first discuss the nature of the fermion sign problem. We consider a fermionic 
path integral 

Z / = ^SignHexp(-5M) (1.1) 

N 

over configurations of occupation numbers n with a Boltzmann weight whose mag- 
nitude is exp(— S[n}) and sign is Sign[n] = ±1. Here S[n] is the action of a corre- 
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sponding bosonic model with partition function 

Z 6 = 5>xp(-S[n]). (1.2) 

[n] 

A fermionic observable 0[n] is obtained in a simulation of the bosonic ensemble as 
(0)f = ^EONSign[n]exp(-5N) = (1.3) 

Z f [n] \ Sl S n / 

The average sign in the simulated bosonic ensemble is given by 

(Sign) = ±-J2 SignM exp(-5[n]) = ^ = exp(-/JVA/). (1.4) 

The expectation value of the sign is exponentially small in both the volume V and 
the inverse temperature f3 because the difference between the free energy densities 
Af = ff — fb of the fermionic and bosonic systems is always positive. Hence, the 
fermionic expectation value (O) f — although of order one — is measured as the ratio 
of two exponentially small signals which are very hard to extract from the statistical 
noise. This is the origin of the sign problem. We can estimate the statistical error 
of the average sign in an ideal simulation of the bosonic ensemble which generates 
iV completely uncorrelated configurations as 

^Sign = \/(Sign 2 ) ~ (Sign) 2 _ exyjpVAf) n ^ 
(Sign) " v^V/Sign) 

The last equality results from taking the large f3V limit and using Sign 2 = 1. In 
order to determine the average sign with sufficient accuracy one needs on the order 
of N = exp(2(3VAf) measurements. For large volumes and small temperatures this 
is impossible in practice. 

In the meron-cluster approach one solves the sign problem in two steps. In the 
first step one analytically re-writes the partition function using new variables such 
that it is possible to exactly cancel all negative weight configurations with configu- 
rations of positive weights. This group of configurations then does not contribute 
to the partition function. The remaining configurations are guaranteed to be pos- 
itive. Thus, in the new variables, one effectively obtains Boltzmann weights with 
Sign = 0, 1. At this stage, despite the fact that all negative signs have been elim- 
inated, only one half of the sign problem has been solved since an algorithm that 
naively generates configurations with Sign = or 1 generates an exponentially large 
number of zero weight configurations. Thus, it is important to introduce a second 
step in which one avoids configurations that have been canceled analytically. In 
practice it is often useful to allow these zero-weight configurations to occur in a 
controlled manner since these configurations may contribute to observables. In a 
numerical algorithm a local Metropolis decision ensures that contributions of and 
1 occur with similar probabilities. This solves the sign problem completely. 
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1.2 The Meron-Cluster Method 



The central idea behind the meron-cluster method is to express the partition func- 
tion, which is originally written as a sum over weights of fermion occupation number 
configurations, as a sum over weights of configurations which contain both fermions 
and new bond variables. Mathematically this means 

Zf = Sign[n] exp(— S[n\) = Sign[n, b] exp(— S[n, b]) (1.6) 

M [n,b] 

where the bond variables b carry information about whether two lattice sites are 
connected or not. A cluster is defined as a set of connected lattice sites whose 
flip exchanges occupied and empty sites, i.e., the occupied sites on that cluster 
are emptied and the originally empty sites now become occupied. This step of re- 
writing the partition function is well known and has been used in earlier attempts 
of constructing fermion cluster algorithms as discussed in and |J. The recent 
progress results from the observation that the sign problem contained in a partition 



function of the form ( |1.6| ) is completely solvable if the magnitude W[n, b] and the 



sign Sign[n, b] of the Boltzmann weight satisfy three requirements: 



1. The magnitude of the Boltzmann weight exp(— S[n, b}) does not change when 
any cluster is flipped. 

2. The effect of a cluster-flip on the sign of the Boltzmann weight Sign[n,6] is 
independent of the orientation of all other clusters. 

3. Starting from any configuration [n,b], it must be possible to flip clusters and 
reach a reference configuration [n re f, b], which is guaranteed to have a positive 
Boltzmann weight. 



In addition, a formula for the change in the sign of a configuration due to a cluster- 
flip has been derived recently [[HJ. Using this formula and other tricks it is usually 
possible to satisfy the first two properties for any given model. Restrictions in the 
class of solvable models arise due to the inability to satisfy the third property listed 
above, i.e., the existence of a reference configuration. However, we have been able 
to construct useful reference configurations for a variety of models. 

The essential consequence of the three basic properties necessary for the meron- 
cluster approach to work is that clusters can be characterized by their effect on the 
sign. Clusters whose flip changes the sign are referred to as merons. Clearly a config- 
uration with a meron-cluster contributes zero to the partition function after a partial 
re-summation over cluster-flips is performed^. All configurations without any merons 
contribute positively to the partition function. Thus the partition function can be 

lr This partial re-summation over cluster-flips is often referred to as an improved estimator. 
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re-written in terms of positive semi-definite Boltzmann weights, i.e., Sign = 0, 1. In 
order to avoid the zero-weight configurations, a re-weighting method that eliminates 
configurations containing multiple meron-clusters is included. In order to perform 
this re-weighting, all observables must also be measured using improved estimators. 
Fortunately, this is possible for most quantities of physical interest. Usually, physi- 
cally interesting observables receive non-zero contributions also from configurations 
containing merons. For example, vacuum condensates receive contributions from 
the one-meron sector whereas observables derived from certain two-point correla- 
tion functions receive contributions from zero- and two-meron configurations only. 
Four-point correlation functions also receive non-zero contributions from four-meron 
configurations. Fortunately, one is usually not interested in multi-point correlation 
functions and for most purposes one can completely eliminate the configurations 
with more than two merons. This is exactly what happens in the re-weighting step 
of the meron-cluster algorithm. In this paper we do not discuss the re-weighting 
step although this is an essential part of the complete solution of the sign problem. 
We assume that once all negative signs are eliminated, one can design a re-weighting 
step that does not generate exponentially long auto-correlation times. 

Although the three requirements for solving the sign problem with the meron- 
cluster method are somewhat restrictive, a large class of Hamiltonians for strongly 
interacting fermions is consistent with them. For example, in the case of spinless 
fermions it is necessary to include a certain amount of nearest-neighbor fermion 
repulsion [Q or a large chemical potential || depending on the hopping term. As we 
will see, similar restrictions arise for fermions with spin. In particular, the standard 
Hubbard model Hamiltonian cannot be treated with our method, while slightly 
modified models in the Hubbard model family can. In this paper we construct a 
large class of SU(2) spin symmetric Hamiltonians for spin one-half fermions with 
nearest-neighbor interactions to which the meron-cluster method can be applied. 
The numerical results for one such model related to s-wave superconductivity have 
been discussed elsewhere HTT |. 



The meron-cluster idea was first developed for a bosonic model with a complex 
action — the 2-d 0(3) model at non-zero vacuum-angle 9 \Y1\. The term meron 



was introduced first to denote half-instantons. Indeed the meron-clusters in the 2-d 



0(3) model are half-instantons | 12| . The instanton number is the integer- valued 
topological charge Q G Z, which gives rise to a complex phase exp(i9Q) that con- 
tributes to the Boltzmann factor. At 9 = n the phase factor reduces to a sign 
exp(z7rOJ = (—1)^, which distinguishes between sectors with even and odd topolog- 
ical charges. In the 2-d 0(3) model the meron-clusters carry half-integer topological 
charges. When such a cluster is flipped, its topological charge changes sign and 
hence the total topological charge of the configuration changes by an odd integer. 
Thus, the phase factor (— 1)^ changes sign when a meron-cluster is flipped. Sim- 
ilarly, a meron-flip in the fermion cluster algorithm changes the sign of a fermion 
world-line configuration. Interestingly a Z(2) topological number can be assigned 
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to each fermionic configuration. Any configuration with Sign = 1 is topologically 
trivial because it represents an even permutation of fermions. A configuration with 
Sign = — 1 is topologically distinct, because an odd permutation cannot be changed 
into an even one by a continuous deformation of fermion world-lines. In this sense, 
configurations with Sign = — 1 are instantons. Flipping any meron changes the 
fermion sign and hence unwinds the instanton. In two spatial dimensions one can 
extend the analogy to the 2-d 0(3) model even further. Then particles can exist 
with any statistics, characterized by a parameter 9 that varies between 6 = for 
bosons and 9 = tt for fermions. For anyons the fermion sign is replaced by a complex 
phase factor exp(i6H), where H e Z is the integer- valued Hopf number. In that 
case the merons are half-Hopf-instantons. 

The merons are not just of algorithmic interest. In fact, they are physical objects 
that allow us to understand the topology of the fermion world-lines and hence the 
origin of the fermion sign. The merons effectively bosonize a fermionic theory. When 
a fermionic model is formulated in terms of bosonic occupation numbers which inter- 
act locally, the fermion sign still induces non-local topological interactions. Hence, 
the resulting theory is not really bosonized. The merons, however, decompose the 
global topology of the fermion world-lines into manageable contributions that are 
local on the scale of the correlation length. This allows us to simulate fermions 
almost as efficiently as bosons. 

1.3 Organization of the Paper 

The paper is organized as follows. In section 2 the essential ideas behind the meron- 
cluster approach are illustrated for the simple case of spinless fermions. We discuss 
how the three properties necessary for the complete solution to the sign problem, 
discussed in section [1.2| , place restrictions on the class of models that are solvable 
with the meron-cluster algorithm. In section 3 we build a class of models of fermions 
with spin which are solvable using meron-cluster methods. We also discuss improved 
estimators for some physical quantities useful for the study of superconductivity that 
receive contributions from zero- and two-meron configurations and present some 
data from a recent numerical simulation. We also show how one can extend these 
models by adding new terms to the Hamiltonian that appear to violate the properties 
necessary for the meron-cluster method to work but do not introduce sign problems. 
In particular, we show how a repulsive Hubbard model becomes solvable in the 
presence of a chemical potential. Section 4 contains a discussion of the efficiency of 
meron-cluster algorithms. Finally, in section 5 we end with some conclusions. 
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2 Meron- Cluster Solution for Spinless Fermions 



Let us illustrate the essential steps that yield solutions to sign problems in the 
meron-cluster approach using the example of spinless fermions. This is essentially 
a combination of results first discussed in 0, || and flQ| . This will serve as a 
preparation for the main subject of this paper — namely fermions with spin. 



2.1 Fermion Path Integral in Fock State Basis 

We consider fermions hopping on a d-dimensional cubic lattice of volume V = L d 
sites with periodic or anti-periodic spatial boundary conditions. The fermions are 
described by creation and annihilation operators cj. and c x with standard anti- 
commutation relations 

{4> C U = °> { c *> c f } = °> (4> c y} = <V ( 21 ) 
Following Jordan and Wigner |T3| we represent these operators by Pauli matrices 

4 = c x = alal-.af^af, n x = c\c x = ~(a? + 1), (2.2) 

where 

of = ± K <] = ^ lm e ijk af. (2.3) 

Here I labels the lattice point x. The Jordan- Wigner representation requires an 
ordering of the lattice points. For example, one can label the point x — (xi,x 2 , ■■■,Xd) 
by I = x\ + x 2 L + ... + XdL d ~ l . It can be shown that the physics is completely 
independent of the ordering. Further, the Jordan- Wigner representation works in 
any dimension. 

Instead of considering the most general model that can be solved using the 
meron- concept, we restrict ourselves to a class of Hamiltonians that conserve particle 
number. From the discussion below it will become clear how to extend the ideas 
presented here to other Hamiltonians while maintaining the ability to solve the 
associated sign problem. With this in mind we focus on a Hamilton operator of the 
form 

H = Y,h x , h (2.4) 

x,i 

which is a sum of nearest-neighbor Hamiltonians h Xj i given by 

Ki = -\{<i c x+i + c l+fx) + U{n x - ^)(n x+t - -). (2.5) 

which couples the fermion operators at the lattice sites x and x + i, where % is a 
unit-vector in the i-direction. Next we decompose the Hamilton operator into 2d 
terms 

H = H 1 + H 2 + ... + H 2d , (2.6) 
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with 



H 



E h 



-i xi H i+d 



E K 



(2.7) 



Xieven 



x=(x 1 ,x 2 ,---,x d ) 
Xiodd 



Note that the individual contributions to a given Hi commute with each other, but 
two different Hi do not commute. Using the Trotter-Suzuki formula we express the 
fermionic grand canonical partition function as 



Z f = Tr{exp[-P(H-^N)]} 

= KmTr {expf-e^ - ^N)\ exp[-e(tf 2 - ^N)}... ex V [-e(H 2d 



2d 



N)} 



M 



(2- 



Here N = J2 X n x is the particle number operator and /x is the chemical potential. We 
have introduced M Euclidean time-slices with e = (3/M being the lattice spacing 
in the Euclidean time-direction. We insert complete sets of fermion Fock states 
between the factors exp[— e(ifj — £N)\. Each site is either empty or occupied, 
i.e. n x has eigenvalue or 1. In the Jordan- Wigner representation this corresponds 
to eigenstates |0) and |1) of of with of |0) = — 10) and crf\l) = |1). The transfer 
matrix is a product of factors 



exp[-e(h X)i - + n x + i))\ = ex P[ £ ( J + ^)] 



x 



/ ex P [-e(f + A)] 



V o 





cosh(f) Ssinh(f) 

Ssinh(f) cosh(f) 











(2.9) 



exp[-c(? - M J 



which is a 4 x 4 matrix in the Fock space basis |00), |01), 1 10) and |11) of two 
sites x and x + %. Here £ = of +1 of +2 ...o^ t _ 1 is a string of Pauli matrices running 
over consecutive labels between / and m, where I is the smaller of the labels for the 
lattice points x and x + % and m is the larger label. Note that £ is diagonal in the 
occupation number basis. 



The partition function is now expressed as a path integral 

Z f = E Sign[n] exp(-S[n]), 



(2.10) 



over configurations of occupation numbers n(x, t) = 0, 1 on a (d + l)-dimensional 
space-time lattice of points (x,t). The magnitude of the Boltzmann factor, 

exp(— S[n\) = 

M-l d 

Y\ n }] exp{-s[n(x,t),n(x + i,t),n(x,t + l),n(x + i,t + 1)]} 

p=0 i=l x=(xi,X2,---,x<i) 

Xieven, t=2dp+i— 1 

x [f exp{-s[n(x,t),n(x + i,t),n(x,t + l),n(x + i,t + l)]}, (2.11) 

Xiodd, t=2dp+d+i— 1 
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Figure 1: A typical fermion world-line configuration in one spatial dimension with 
periodic boundary conditions. The product of £ factors from each plaquette is nega- 
tive for this configuration since the two fermions exchange their positions with each 
other. 

is a product of space-time plaquette contributions with 

exp(-s[0, 0, 0, 0]) = exp[-e(| + ^)], 

exp(-s[0, 1, 0, 1]) = exp(-s[l, 0, 1, 0]) = cosh(|), 

exp(-s[0, 1, 1, 0]) = exp(-s[l, 0, 0, 1]) = sinh(|), 

exp(-s[l, 1, 1, 1]) = exp[-e(| - A)]. ( 2 .12) 

All the other Boltzmann factors are zero, which implies several constraints on allowed 
configurations. Note that we have dropped the trivial overall factor exp[e( ¥ + £)] 
in eq. (|2.9| ). The sign of the Boltzmann factor Sign[n] also is a product of space-time 
plaquette contributions sign[n(x, t), n(x + i, t),n(x, t + 1), n(x + i, t + 1)] with 

sign[0, 0, 0, 0] = sign[0, 1, 0, 1] = sign[l, 0, 1, 0] = sign[l, 1, 1, 1] = 1, 
sign[0,l,l,0] = sign[l,0,0,l] = S. (2.13) 

It should be noted that the non-local string of Pauli matrices S gets contributions 
from all lattice points with labels between I and m. This would make an evaluation 
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of the fermion sign rather tedious. Also, it is not a priori obvious that Sign[n] is 
independent of the arbitrarily chosen order of the lattice points. Fortunately, there is 
a simple way to compute Sign[n] which is directly related to the Pauli principle and 
which is manifestly order-independent. In fact, Sign[n] has a topological meaning. 
The occupied lattice sites define fermion world-lines which are closed around the 
Euclidean time-direction. Of course, during their Euclidean time-evolution fermions 
can interchange their positions, and the fermion world-lines define a permutation 
of particles. The Pauli principle dictates that the fermion sign is just the sign of 
that permutation. When one works with anti-periodic spatial boundary conditions, 
Sign[n] receives an extra minus-sign for every fermion world-line that crosses the 
boundary. Figure p] shows a typical configuration in one spatial dimension along 
with examples of contributions to the Boltzmann weight arising from local space- 
time plaquettes. 



2.2 Observables in Fermionic Fock Space 

A variety of observables of physical interest take on a simple form in Fock space. 
Typically the expectation value of an operator O is given by 

(°)f = ^E°MSign[n]exp(-S[n]). (2.14) 

Z f [n] 

For example, the total particle number relative to half-filling 

# M - j = E ("OM) " 5) , ( 2 - 15 ) 

which is conserved and hence the same in each time-slice, and its disconnected 
susceptibility 

*»-£(("-?)'), < 2 - 16 ' 

are both operators that are diagonal in the occupation number basis. Other quan- 
tities of physical interest are the staggered occupation 

[n] = e J2(-ir +X2+ - +Xd (n(x,t) - I) , (2.17) 



and the corresponding susceptibility 



Xo = ± ((O 2 ), - (O)}) . (2.18) 
We can also measure the winding number Wi[n] around the spatial direction i 

M-l ^ 
Wi[n\ = ^2 E ~Jj {fi[n(x,t),n(x+i,t),n(x,t+l),n(x+i,t+l)], [1,0,0,1] 

Xieven, t=2dp+i— 1 
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3{n(x,t),n(x+i,t),n(x,t+l),n(x+~i,t+l)],[0,l,l,0} } 



+ E 7 {d{n(x,t),n( 



\ u [n(x,t),n(x+i,t),n(x,t+l),n(x+i,t+l)],[l,0,0,l] 
Xjodd, t=2dp+d+i-l 



u [n(x,t),n(x+i,t),n(x,t+l),n(x+i,t+l)],[0,l,l,0] 

}. (2.19) 



This counts 1 for a plaquette configuration [1,0,0,1], i.e., for a fermion hopping 
in the positive i-direction, and —1 for [0, 1, 1,0], i.e., for a fermion hopping in the 
opposite direction. As a result, Wi[n] counts how many fermions wrap around the 
i-direction during their Euclidean time-evolution. Again, one can define a corre- 
sponding susceptibility 

XWi = P(Wi)f- (2-20) 

We will use a similar quantity as an order parameter for superconductivity in models 
of fermions with spin. It measures the response of the system to a twist in the 
spatial boundary conditions. In the infinite volume limit the system feels the spatial 
boundary only in a phase with long-range correlations. 

Observables like the two-point Green function defined by 

Tr {exp [-(/? - t)H] c.exp [-(* - t')H] ctexp [-t'H]} 
^Wj- Tr {exp [-/?#]} { ' 

are quite useful to extract the spectral information of H . These expectation values 
of non-diagonal operators can also be obtained like those of the diagonal operators 
using 

G(x,t;y,t') = -L£Sign[n']exp(-£K]). (2.22) 

Z f [n'} 

where the configurations [n'\ that contribute are different from the configurations 
that contribute to the partition function due to the violation of fermion number 
at the space-time sites (x,t) and (y,f). In order to allow for these violations it is 
convenient to introduce two new time-slices at t and t'. The factors exp(—S[n']) and 
Sign[n'] can be calculated in the same way as before on all the time-slices except on 
the two slices where fermionic creation and annihilation operators are introduced. 
At those two slices the fermion occupation number changes appropriately and the 
Jordan- Wigner representation can be used to figure out the extra sign factor coming 
from the string of 03 associated with fermionic creation and annihilation operators. 



Interestingly, as discussed in [14]], once the path integral is re- written in terms of 
cluster variables, it is straightforward to keep track of these new configurations 
along with the configurations that contribute to the path integral and hence such 
observables can also be measured. 
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2.3 Cluster Decomposition of the Partition Function 



Up to now we have derived a path integral representation for the fermion system in 
terms of bosonic occupation numbers and a fermion sign factor that encodes Fermi 
statistics. The system without the sign factor is bosonic and is characterized by the 
positive Boltzmann factor exp(—S[n}). Here the bosonic model turns out to be a 
quantum spin system with the Hamiltonian 

H = YX^SlSU + + US 3 x S 3 x J -fiJ2 S l (2-23) 

x,i x 

where S x = \o\ is a spin 1/2 operator associated with the lattice site x that was la- 
beled by /. The case U = t corresponds to the isotropic anti-ferromagnetic quantum 
Heisenberg model, while U = represents the quantum XY-model. The chemical 
potential plays the role of an external magnetic field. 

It is well known that the partition function for the above quantum spin model can 
be re-written in terms of cluster variables. This cluster representation is at the heart 
of the the loop-cluster algorithm |T5|, [TB[ that is used to simulate the quantum spin 
model very efficiently. Interestingly, the meron-cluster algorithm for the fermionic 
model is also based on a similar cluster decomposition of the partition function. 
The essential differences come from the fermion sign factor and are encoded in the 
topology of the clusters. Once the connection between the fermion permutation 
sign and the topology of the clusters is understood, the cluster algorithm for the 
quantum spin system can be modified easily to deal with the fermionic model. 

Let us first discuss the essential steps in the cluster decomposition of the partition 
function of the quantum spin model. The main idea is to take a configuration of 
quantum spins and re-write the weight as a sum over weights of configurations with 
spins and bonds, such that the magnitude of the weight of the configuration does 
not change when the set of spins connected by bonds are flipped. For the problem 
at hand this decomposition can be accomplished at every plaquette. For example, 
the Boltzmann weights of each plaquette given in eq. ( |2.12| ) can be written as 

exp(— s[n(x, t),n(x + i, t),n(x, t + 1), n(x + i, t + 1)]) = 

exp (— s[n(x, t), n(x + i, t),n(x, t + 1), n(x + i, t + 1); b[\ , (2.24) 

b=A,B,C,D 

where b = A,B,C,D represent four types of bond break-ups. The Boltzmann 
weights exp (— s[n;b}) can be arbitrary as long as they satisfy eq. ( |2.24| ). Here we 
choose 

exp (-s[0, 0,0,0; A}) = exp (— s[l, 0, 1, 0; A}) = exp (— s[0, 1, 0, 1; A}) = 

exp (— s[l, 1, 1, 1; A]) = W A = ^ jexp (-e [| + +exp (-|)} , 
exp (-s[0, 0,0,0; B}) = exp (-s[l, 0, 0, 1; B]) = exp (-s[0, 1, 1, 0; B]) = 
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-exp B]) = W B = _exp(~) sinh , 

exp (-s[0, 0, 0, 0; C]) = exp (s[l, 0, 0, 1; C}) = exp (-s[0, 1, 1, 0; G\) 

exp(-s[l,l,l,l;C]) = W c = \ {exp (-e [| - ^]) - exp (- 
exp(-s[l,0,l,0;L>]) = exp (-s[0, 1, 0, 1; £>]) = exp (— s[0, 1, 1, 0; D]) 

exp (-s[l, 0,0, !;£>]) = W D = X - {exp (|) - exp 



h — 

2 2d 



^2.25) 



to be the only non-zero weights. The above equations are shown pictorially in 
figure pi. If we interpret the spin states as fermion occupation numbers, clearly 



W[n, b] satisfies the first property discussed in section [L^. We will now assume 
that [A < and t + fi/d > U, so that all the Boltzmann weights exp(— s[n; b}) of 
spin and bond configurations are positive. However, notice that there is an extra 
negative sign for each plaquette with configuration [1, 1, 1, 1; B\. Such negative signs 
must be considered separately as a contribution to the sign of a given spin and 
bond configuration. Explicitly this is given by (— 1) Nb where Nb is the number 
of plaquettes of the type [1, 1, 1, 1; B]. Generally, such extra sign factors must be 
avoided since they can lead to sign problems. However, since we anyway have to 
deal with the non-local fermion signs, which we have suppressed until now, we will 
allow for such extra negative sign factors. As we will see below, these local sign 
factors can sometimes be helpful in canceling sign factors that arise due to fermion 
permutations. 



Based on the decomposition given in eq.( 2.24|) the original partition function of 
the fermionic model can be written as 

Z f = J2 Sign[n,6] exp(-S[n, 6]), (2.26) 

[n,6] 

a sum over configurations of occupation numbers n(x, t) — 0, 1 and bond variables 
b = A, B,C, D on each plaquette. The sign Sign[n, b] is the product of the original 
fermion permutation sign due to the occupation numbers and (— 1) Nb . The new 
configurations naturally connect neighboring sites through the bonds that live on 
each individual space-time interaction plaquette as represented pictorially in figure 
A sequence of connected sites defines a cluster. Since each bond on a plaquette 
connects two sites, every cluster is a simple closed loop. Thus we have written 
the partition function of a fermionic model in terms of the statistical mechanics of 
closed interacting loops with not necessarily positive Boltzmann weights. Although 
we have illustrated this using a particular model, it is straightforward to repeat the 
above steps for any model. When a cluster is flipped, the occupation numbers of 
all the sites on the cluster are changed from n(x,t) to 1 — n(x,t). In other words, 
a cluster-flip interchanges occupied and empty sites. The weights given in eq. (|2.25|) 



guarantee that the magnitude of the weight of a configuration does not change when 
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exp(-s[0,0,0,0]) 



exp(-s[l, 1,1,1]) 



exp(-s[0, 1,1,0]) 



exp(-s[ 1,0,0,1]) 




exp(-s[l,0,l,0]) 




exp(-s[0, 1,0,1]) 
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Figure 2: The pictorial representation of the decomposition of plaquette weights of 
fermion occupation number configurations in terms of weights of bonds and fermion 
occupation numbers. 
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-Wb 



Figure 3: A typical configuration of fermion occupation numbers and bonds in one 
spatial dimension. 

a cluster is flipped. A typical configuration in one spatial dimension is represented 
schematically in figure 0. 



2.4 Meron-Clusters and Reference Configurations 

Let us now consider the effect of a cluster-flip on the sign factor. In the example 
considered, the sign of a given configuration of spins and bonds arises from two 
factors. One is the permutation of fermion world-lines which is independent of the 
bond configuration and the other is due to plaquettes of type [1, 1, 1, 1; B] each of 
which contributes a negative sign to the Boltzmann weight. Obviously, a cluster-flip 
either changes the sign of a configuration or not. In general, the effect of the flip 
of a specific cluster on the fermion sign depends on the orientation of the other 
clusters. For example, a cluster whose flip does not change the sign now, may very 
well change the sign after other clusters have been flipped. In other words, in general 
the clusters affect each other in their effect on the configuration's sign. Recently, 
a formula for how cluster-flips affect the fermion permutation sign was discovered 
[10|1 . The relevant information is encoded by an integer which can be determined 



by starting at an arbitrary point on the loop-cluster and traversing around it in 
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either direction and nipping the sites as they are encountered. As we go around the 
loop, we denote the number of horizontal bonds (D-type bonds) encountered while 
entering an occupied site from an empty site (before they are flipped) by N^, the 
number of cross-bonds (B- or C-type bonds) traversed while the other cross-bond 
on the same plaquette connects two occupied sites by N x , and the temporal winding 
of the cluster by N w . Then the cluster-flip changes the fermion permutation sign if 
and only if 

N = N h + N w + N x (2.27) 



is even. For details of how this can be derived we refer the reader to ]TD[ . A similar 
formula can be derived for the signs associated with the plaquettes of the type 
[1, 1, 1, 1; B]. In this case, again as we traverse the loop-cluster and flip the sites, if 
N x b is the number of -B-type cross-bonds encountered when the other cross-bond 
connects two occupied sites, then the cluster-flip changes the sign associated with 
the plaquettes of the type [1, 1, 1, 1; B] if N x b is odd. Using this formula we discover 
that in the present model the effect of a cluster-flip on the sign of a configuration 
in general depends on the orientation of other clusters. In particular, if two clusters 
cross each other in an odd number of C-type cross-bonds then the flip of one cluster 
has an effect on whether the other cluster-flip changes the configuration sign or 
not. Interestingly, once we forbid C-type cluster break-ups, i.e. when Wc = 0, the 
clusters have a remarkable property with far reaching consequences: each cluster 
can be characterized by its effect on the fermion sign independent of the orientation 
of all the other clusters. We refer to clusters whose flip changes Sign[n, b] as merons, 
while clusters whose flip leaves Sign[n, b] unchanged are called non- merons. The flip 
of a meron-cluster permutes the fermions and changes the topology of the fermion 
world-lines. Notice that the -B-type cross-bonds are still allowed thanks to the extra 
negative signs that arise due to signs from [1, 1, 1, 1; B] plaquettes. Thus we have 
been able to ensure the first two properties discussed in section |1.2| , provided we 
satisfy the constraints [i<0,t + fi/d>U and Wc = which implies t = U — \xj 1 d. 

When the cluster-flips affect the sign factor independently, it is easy to average 
analytically over all the cluster-flips. Such an average over cluster-flips for a given 
bond configuration is referred to as an improved estimator in the language of cluster 
algorithms. We will discuss improved estimators of some observables in the next 
section. Averaging the sign over cluster-flips, the fermionic partition function can 
be re-written as 

Z S = Signl6] exp(-S[n, &]), (2.28) 

[n,b] 

where 

SPI = ^ E Sign[n,6]. (2.29) 

" cluster— flips 



Thus Sign [b] = if at least one of the clusters is a meron. In order to solve the 
sign problem it is important to eliminate configurations with meron- clusters from 
the Monte Carlo sampling. Unfortunately, this alone does not guarantee that the 
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contributions from the zero-meron sector will always be positive, although it may 
alleviate the sign problem considerably. With no merons in the configuration it is 
clear that the sign of a configuration remains unchanged under cluster-flips, but one 
could still have Sign[n,6] = — 1. In order to solve the sign problem completely it 
is necessary to show that configurations with no merons are always positive. This 
takes us to the third property discussed in section fL2| , the need for a reference 
configuration with a positive weight. Interestingly, there exists a class of models 
where this is achievable. It is then easy to show that Sign [b] = 1 in the non- 
meron sector. This solves the sign problem completely. Thus the fermionic partition 
function is obtained as a sum over weights of configurations in the zero-meron sector, 



where W[b] is obtained as a product of local weights Wa, Wb and Wd- 

There are at least two types of reference configurations that arise for spinless 
fermions. For example, in addition to setting Wc = 0, if we eliminate the 5-type 
bonds completely by also setting Wb = 0, the remaining cluster break-ups have a 
remarkable property. They guarantee that sites inside a cluster obey a pattern of 
staggered occupation, i.e. either the even sites (with X\ + %2 + ... + Xd even) within 
the cluster are all occupied and the odd sites are all empty or vice versa. This 
guarantees that the clusters can be flipped such that one reaches the totally staggered 
configuration in which all even sites at every time-slice are occupied and all odd sites 
are empty. In this half-filled configuration all fermions are static, no fermions are 
permuted during the Euclidean time-evolution. Further, since the model has no 
other sign factors, we know that Sign[n re f , b] — 1. Since any configuration can be 
turned into the totally staggered configuration by appropriate cluster-flips, property 
three is indeed satisfied. The condition Wb = and Wc = leads to the restriction 
fj, = and t = U . This leads to an interesting interacting spinless fermion model. A 
variant of this model has been studied in 0, [5|, |6|, 0] where it has been shown that 
the model can be used to study the spontaneous breaking of a Z(2) symmetry. 

The second reference configuration is obtained by setting Wd = along with 
Wc = 0. In this case we see that all sites in a given cluster have the same occupation 
number, i.e., all are either empty or all are occupied. This means that all clusters can 
be flipped such that one reaches a configuration where all sites are empty. Since there 
are no fermions at all, there is no fermion permutation sign in this configuration. 
Further, there are no negative signs because of plaquettes of the form [1, 1, 1, 1; B], 
since all sites are empty and [0, 0, 0, 0; B] is positive. Therefore the configurations 
with no merons necessarily have a positive sign. The condition Wc = and Wd = 
leads to the restriction U — and t = —fi/d with [i < 0. This model leads to free 
non-relativistic fermions and has been discussed in ||. 

Figure f| shows a typical configuration and the associated reference configuration 
for each of the two classes discussed above. The reference configurations are closely 




£ 2"'W[b], 
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Figure 4: Figure (a) shows a typical fermion and bond configuration without cross- 
bonds and (b) shows the corresponding reference configuration which has a positive 
sign. Figures (c) and (d) show similar configurations for the case without horizontal 
bonds. 



connected with the ground state properties of each model. The staggered reference 
configuration naturally describes theories with spontaneously broken ground states 
whereas the totally empty configuration leads to free non-relativistic fermions. It is 
possible to include a limited form of additional interactions without destroying any of 
the useful properties in these two classes of theories. However, it would be interesting 
to find new reference configurations which can lead to interesting physics. Of course, 
there are a variety of non-translationally invariant reference configurations which 
could describe electronic properties of certain crystalline structures. These models 
may be worth exploring. As we will show below, novel reference configurations also 
arise when spin is introduced. 
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2.5 Building-Blocks of Cluster Models 



The results of the last two sections can be used to synthesize simple rules that help 
in building models of spinless fermions which automatically satisfy the first two 
properties of the meron-cluster approach, i.e. the weight of a configuration of clusters 
does not change under a cluster-flip and the effect on the sign of a configuration due 
to a cluster-flip is independent of all other clusters. In the models constructed 
above we started with four types of bond interactions namely the A-, B-, C- and 
D-types as shown in the figure |5|. Each of these interactions can be associated with 
an operator acting on the two-site Hilbert space. For example, a plaquette with 
A-type (vertical) bonds gives the same positive weight to all diagonal elements. 
Hence A-type bonds are equivalent to the unit matrix. Similarly, the B-, C- and 
-D-type bonds are equivalent to operators tabulated in table |I|. The fermionic nature 



of the creation and annihilation operators are encoded in the formula of eq.( [2.27| ) 
which is used to figure out the relative sign of configurations when clusters are 
flipped. Using eq.( [2.27| ) along with the weights of the building-blocks, we found that 





B 



D 



Figure 5: Building-blocks of the two-site interactions discussed in the text. 



Bond-type 


Transfer Matrix 


A 




1 




B 




C x Cy H~ CyC x Ti x fly 


\- 1 


C 


c x c y 


+ dc x + 2(n x - 1/2) (n y - 


1/2) + 1/2 


D 




+ c\c x - 2(n x - 1/2) (n y - 


1/2) + 1/2 



Table 1: The two-site building-blocks of figure^ and their transfer matrix contri- 
butions. 

type C bonds are forbidden as they violate the required property that cluster-flips 
should have an independent effect on the sign. Thus the A-, B- and D-type bonds 
provide the basic building-blocks for constructing Hamiltonians which satisfy the two 
necessary requirements of the meron-cluster approach. We only need to ensure the 
existence of a reference configuration. This restriction then forbids models which 
simultaneously allow B- and D-type bonds. Thus, we find two types of solvable 
models discussed above. It is easy to construct the Hamiltonian of these models 
using the corresponding building blocks. For example, the model with A- and 5-type 
bonds leads to the two-site transfer matrix T xy = Wa + Wb [ c x c y + c y c x — n x — n y + 1] , 
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which in the limit e — > yields T xy = 1 + e[cj.c y + c y c x — n x — n y + 1). Since all 
two-site interactions contribute equally to the full Hamiltonian, H = J2< xy > h xy , in 
this limit we can identify T = 1 — eH, which yields 

h xy = -{clcy + CyC x ) + n x + n y ( 2 -31) 

up to additive and multiplicative constants. Similarly the model obtained with A- 
and -D-type bonds leads to 

h xy = -{clc y + c\c x ) + 2{n x - l/2){n y - 1/2). (2.32) 

Thus, one can construct the Hamilton operator once the building-blocks are chosen 
in terms of the bond break-ups on an elementary plaquette. Later, we will use this 
procedure to construct models of fermions with spin. In eq. ( |2.31| ) and eq. fl2.32| ) 



we have chosen to normalize the hopping term and set t = 1. In the remainder of 
this article we will continue to adopt this normalization for convinence. 



2.6 Improved Estimators 

Based on the previous discussion it is clear that the partition function can be written 
entirely in terms of cluster variables as 

z f = J2 2N ° w \ b \ ( 2 - 33 ) 

[b],zero— meron 

where Nc is the total number of clusters and the sum is over cluster configurations 
without meron-clusters. Similarly, it is possible to find expressions for observables 
in the language of clusters. For example, consider the deviation of the average 
occupation number of a configuration from half-filling which is defined as 

An = n -\ = v^{ nx ~l)- (2 - 34) 

Clearly it must be possible to sum over x belonging to a particular cluster C e [b] 
first and then sum over all possible clusters. If we define 

= X>* - 1/2), (2.35) 

it is easy to see that An = (l/V)J2c^c- Performing a partial average over flips 
of fermionic occupation numbers belonging to various clusters one finds that this 
average is non-zero only in the one-meron sector and in this sector 

( An )ciustcr-flip S = 77 ^e mcro n (2.36) 
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where fic mcron is to be calculated for the meron-cluster which is flipped into the 
reference configuration. Thus the complete answer turns out to be 

(An) = — J2 ^f^2 Nc W[b\. (2.37) 

f [ft], one— meron 

In the case of spinless fermion models discussed above, it is easy to see that flc is 
equal to the temporal winding of the cluster C up to a factor of 2. For the model 
described by the Hamiltonian of eq.( 2.32 ), (An) = since there is no contribution 
from the single-meron sector to the partition function. Similar improved estimators 
for the chiral condensate || and susceptibility || can be constructed. In section |3l^ 
we will consider improved estimators for the susceptibilities relevant for the study 
of superconductivity. 



3 Models of Fermions with Spin 

There are many ways to construct models of fermions with spin such that the three 
properties necessary for the meron-cluster method to work, as discussed in section 
|1.2| , are satisfied. The most general approach is to start from the Hamilton operator 
of interest and proceed with the same steps as discussed in the earlier section. 
If the Hamilton operator is a sum of only two-fermion interactions, then nothing 
qualitatively new emerges since it is possible to consider the spin degree of freedom 
as two "spatial" layers representing up and down spins. All on-site interactions that 
can occur between the two spin layers can be introduced on a separate time-slice. 
However, it is now possible to find new reference configurations which are physically 
meaningful and new models can be constructed using them. 

On the other hand, it is natural to allow four- "site" interactions which connect 
nearest spatial neighbors and the two spin layers. In this case new situations not 
considered in the previous sections can arise. As an example, consider the transfer 
matrix element given in figure |^. Although such four-site interactions have not been 
considered earlier, they can be thought of as two correlated two-site interactions. 
For example, the shown interaction can arise either due to a "correlated" fermion 
hop which conserves spin or a "correlated" spin-flip. These two alternatives are 
shown on the right-hand side in the same figure. If we assume that these new 
interactions arise from such correlated two-site interactions, then it is possible to 
use the technology developed for spinless fermions. Here we will take this approach 
and consider models with four-site interactions. 

Instead of trying to find the most general solvable model, we again restrict our- 
selves to a class of models that can be built with simple correlated two-site interac- 
tions. In the next section we will enumerate some of the simplest correlated two-site 
interactions and in the following section we will build one of the simplest solvable 
models which has the symmetries of the Hubbard model. 
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Spin down Fermions 




Figure 6: Configurations with four-site interactions can arise as correlated two- 
site interactions in different ways. The two possible interactions result in different 
fermion permutation signs. 
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Figure 7: Spin- conserving bond break-ups that form the building-blocks of models 
discussed in the text. 

3.1 Building-Blocks of a Cluster Model 



Let us first discuss the simplest building-blocks for constructing models of fermions 
with spin that can be solved using meron-cluster techniques. In the case of spinless 
fermions, these building-blocks were discussed in section |2.5| . There it was shown 
how these elementary blocks satisfy two of the necessary properties and how one 
can build models with a reference configuration that ensures that all the criteria 
of solvability are met. We also discussed how it is possible to find the Hamilton 
operator for a model that is solvable using the meron-cluster approach. In the case of 
fermions with spin the various bond break-ups (E, F, G, H, I, J) which conserve spin 
are shown pictorially in figure [7| In table ^| the transfer matrix elements associated 
with each of the building-blocks is tabulated. Using the building-blocks it is possible 
to find models which satisfy all the three required properties of solvability discussed 



in section 1.2. As an illustration we construct a Hubbard- type model below. 
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Bond-type 



Transfer Matrix 



[cxi^Cyi + c yl ^c xl - n xi -n yl + l 



E 

F 
G 
H 

I 

J 



. . .. .. . n vl 

X [Cx^Cyl + Cy\ ] Cx^ - ™:rT - U y i; + 1] 

Cxl'Cyl + c yl ] °x[ - n xi - n yl + 1 
Cx^Cyl + Cy^c x ^ - - n y <t + 1 



[cxi*c yl + cyjci - 2{n xi - 1/2) - 1/2) + 1/2] 
<[c^ t c 2/T + cy^cxt - 2(n* T - l/2)(ny T - 1/2) + 1/2] 
Cx[ ] c yi + c^c^ - 2(n xi - l/2)(n v i ~ V 2 ) + V 2 
Oct^t + c i/T tc ^T ~ 2 ( n ^T ~ V 2 )( n yT ~ V 2 ) + V 2 



Table 2: T/ie operators associated with the building-blocks of spin conserving Hamil- 
tonians shown in figure 0. 

3.2 A Hubbard-Type Model 

The Hubbard model is one of the simplest models used to describe superconductivity. 
In particular, it is believed that the 2-dimensional repulsive Hubbard model may 
describe the physics of high-T c cuprate superconductors flTj . The Hamilton operator 
of the Hubbard model is given by 



H 



<xy>,s=1,l 



n 



x\ 



2 K nxl ~ 2 



/i(n xT + n xi ) 



(3.1) 



This model has a global SU(2) spin symmetry whose generators are given by 

S+ = E c U c *l> S ~ = E c U c *h S 3 = l EKt - n xi) ( 3 - 2 ) 

X X X 

and a global SU (2) pseudo-spin symmetry at \i — whose generators in two dimen- 
sions are given by 



</ + = E(-i 



>xi+a;2„t J 



C 'x\ C 'xii J 



E( 



\x ± +x 2 



Cx\C xh J 3 = - E(?M + U xl ~ K 



(3.3) 

The chemical potential breaks the pseudo-spin symmetry to a £7(1) fermion number 
symmetry. Superconductivity is a result of the spontaneous breaking of this U(l) 
symmetry. In two dimensions, due to the Mermin- Wagner theorem, it is believed 



that this symmetry breaking must follow the Kosterlitz-Thouless predictions [18 



For U < the theory is expected to show s-wave superconductivity at low tempera- 



tures |19], £(J. Then the model can be formulated such that there is no sign problem 
in traditional quantum Monte Carlo algorithms. However, due to the difficulty of 
these conventional fermion algorithms the determination of the critical temperature 
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Spin up Layer 



Spin down Layer 



Figure 8: A typical reference configuration for the Hubbard-type model described in 
the text. 



is still controversial [EIJ. The repulsive model with U > 0, on the other hand, is 
more interesting but suffers from a sign problem away from half-filling 

Variants of the Hubbard model are equally interesting. For example, the t-J 
model has also been extensively studied. Here we refer to such models as Hubbard- 
type models. Let us now construct one such model which is solvable using meron- 
cluster techniques. The model is based on the four-site transfer matrix elements E 
and H shown in figure [?[ In this model the clusters are identical in the spin-up 
and the spin-down layers and any configuration of clusters can be flipped into a 
configuration which is identical in the two layers. Such a configuration always has 
a positive sign and can hence play the role of a reference configuration. A typical 
reference configuration is shown in figure ^. Interestingly, unlike in the spinless 
fermion case, the reference configuration is non-static and not unique. This leads 
to interesting physical consequences. Another feature of the model is that there are 
always an even number of meron-clusters: if the cluster in the spin-up layer is a 
meron then so is the cluster in the spin-down layer. 

Using the operators of table |2| the nearest-neighbor Hamilton operator for the 
model turns out to be 



-[cjc yl + Cyjc xl - 2(n xl - l/2)(n yl - 1/2) + 1/2] 
x [c x ^Cy^ + Cy^c^ - 2(n* T - 1/2) (n,,! - 1/2) + 1/2] 
- A [c x ^c yl + c yi ] c xi - n xl - n yl + 1] 

X [Cx^Cy^ + Cy^CxI - - Uy^ + 1]. 

which leads to the Hamilton operator H = J2< X y> ^xy The operator given in eq. 
can be simplified to Wl 



(3.4) 



h 



■'-y 



■'xs ^ys 



CyS Cy S ^ C X g){n X y \^){n x y 3 ) -|- 2,(S X • Sy ~\- J X • t/y) 
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- 4 [n 



2 J 1 I 



^1 



2) + 2(5, -Sy + J'x- J y ) - 4J, 3 J. 3 



v1 



(3.5) 



"I - A (Ors '--ys -|- Cy S ^ C X s) ij^xy 
- {^-1)^1- l 

where we have used the definition n xy = n x ^ + n x i + n y i + n y i- The operators S x 
and J x are given by the expressions that depend on the site x inside the summation 
signs in eq.( [3.2|) and (|3~ 



The symmetries of the Hamilton operator given in eq. ( |3.5|) can easily be found. 
It is interesting to note that the model is invariant under the SU(2) spin symmetry 
and the U(l) fermion number symmetry. Further for A = the £7(1) symmetry 
enhances to the full SU(2) pseudo-spin symmetry. Although unconventional, the 
physics of this model is closely related to that of the usual attractive Hubbard model. 
In particular, it was found recently that this model shows s-wave superconductivity 
even at half-filling ||11||. Later we will see that it is possible to add a chemical 



potential and an on-site attractive interaction to the above Hamiltonian without 
introducing sign problems. Further, when A = 0, it is possible to include an on-site 
repulsion between the fermions of opposite spin even in the presence of a chemical 
potential as long as a certain inequality is satisfied. It is interesting to investigate if 
this "repulsive" model can be a useful candidate for high-T c superconductivity. 



3.3 Improved Estimators 

Since the model satisfies all the properties of the meron-cluster approach, its parti- 
tion function can be expressed purely in terms of clusters, 

Z = J2 ^ W[b], ( 3 - 6 ) 

[b] ,zero— meron 

where Nq is the number of clusters on each spin layer. Since the clusters are identical 
we get an extra factor of two as compared to eq.( |2.3U| ). This novel way of re- writing 
the partition function involves a partial re-summation over cluster-flips and can 
be extended to other observables as well. We discussed how the deviation of the 
occupation number from half-filling can be expressed in this new language in section 



Let us consider observables relevant for s-wave superconductivity. An important 
observable is the s-wave pair susceptibility 

1 rP 
j3VZ Jo 

withp + = J2x c it c lj. the pair creation andp~ = (p + Y the pair annihilation operators. 
An improved estimator for this susceptibility is easy to construct using the results 



Pl = f Q dt Tr [ exp[-09 - t)H] (p+ + p") exp[-tH] (p+ + p") ] (3.7) 
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of [FJ]. In particular, one can show that non-zero contributions to the two-point 
correlations arise only if the creation and annihilation operators occur on the same 
cluster. Further, in the case of the pair susceptibility it is also important that the 
clusters these operators reside on have been flipped to be identical. The contribution 
is then proportional to the sum of all possible creation and annihilation points along 
the cluster. This just gives a factor of the size of the cluster squared. Averaging 
over the different possible cluster-flips gives for the zero-meron sector 

("^^clustcr— flips, zero— moron N^V ^ ' 2^^ ' (^"^0 

where \C\ is the size of cluster C. In the two-meron sector we get 

("^k) cluster— flips, two— meron jy2 y 2 l^ meron I ' 

(3.9) 

Note that in the two-meron sector there are two identical meron-clusters: one in the 
spin-up and one in the spin-down layer. 



Another observable relevant to superconductivity is the helicity modulus |f23| , |24 
which is defined in terms of the winding number as 

r = ^(w 2 x + w*). (3.10) 

Here W x (W y ) is the total number of fermions winding around the spatial bound- 
ary in the x- (y-) direction. To measure W x we need to examine each interaction 
plaquette crossing the boundary in the x-direction. We count +1 for each fermion 
(up or down spin) that hops across the plaquette in the forward direction, —1 for 
each that hops in reverse and otherwise. This winding for a plaquette can then 
be easily expressed in terms of the occupation numbers of the sites touching the 
plaquette. Consider first just the up spin layer. The winding number on a plaquette 
that originates from the point x, t and goes in the forward 2-direction for one spin 
layer is 

W x ^ = - (n Xtt + n x+lt+1 - n x+u - n x> m ) . (3.11) 

We can also replace each n above with n — 1/2 without changing the results to 
make the values more symmetric. Now every site counts as ±1/4 (on a single spin 
layer) depending on its occupation number and where it is located on the boundary 
plaquette. Each cluster can then be assigned a spatial winding number which is just 
the sum of the weights associated with each point on the cluster. Flipping a cluster 
inverts all the occupation numbers and thus changes the sign of the spatial winding 
number for that cluster. The total winding number in the x-direction is then a sum 
over all clusters and both spin layers 

W X = Y.W C J + W C J . (3.12) 

c 
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After squaring this we average over all cluster-flips. In the zero-meron sector one 
obtains 

(^ 2 ) clustcr - fli ps, z c r o- m c r on = £(W?) a + = 2 £(IUf ) 2 . (3.13) 

c c 

In the two-meron sector we get 

(W x 2 ) cluster _ flips , two _ meron = 2 \Wj^^\Wjr* m L\ = 2(W£—^ 2 . (3.14) 

Again, the two-meron sector consists of one meron-cluster in the spin-up layer and 
one meron-cluster in the spin-down layer both of which are identical. 

Using similar procedures it is also possible to calculate other observables like the 
spin susceptibility which is important to detect magnetic transitions, the d-wave 
pair susceptibility which is important for high-T c superconductivity, and various 
other two-point correlation functions which are of interest. In principle, it is also 
possible to calculate higher-point correlation functions. However, these will become 
increasingly noisy and would require more statistics. 



3.4 Kosterlitz-Thouless Transitions 

As discussed above, the attractive Hubbard model has been studied extensively 
as a toy model for s-wave superconductivity |20, pif. It is expected that below a 



critical temperature transportation of fermion number through the bulk becomes 
easy, leading to superconductivity (or more appropriately superfluidity since the 
symmetry is not gauged in the model). In higher dimensions this is related to the 
spontaneous breaking of the U(l) fermion number symmetry. In two dimensions, 
since this is forbidden due to the Mermin- Wagner theorem, superconductivity occurs 
due to the Kosterlitz-Thouless (K-T) phenomena p 



A variety of striking predictions exist for a K-T transition. For example, at 
temperatures above the superconducting transition temperature T c , the pairing sus- 
ceptibility defined in eq.( |3.7| ) should reach a constant for large enough volumes. As 
T c is approached from above this constant should diverge as exp(a/ y/T — T c ). Below 
T c one should be in a phase with long-range correlations such that the finite- volume 
pair susceptibility diverges as 

P L oc L 2 ~^ T \ (3.15) 

The critical exponent rj equals 1/4 at T c and continuously changes to zero as the 
temperature approaches zero. Unfortunately, none of these have ever been confirmed 
using conventional algorithms even in the attractive Hubbard model which does not 
suffer from a sign problem in the conventional formulation. 



We propose that the model given in eq. ( |3.5| ) is useful for studying the qualitative 
physics of the attractive Hubbard model since it provides a computational advan- 
tage. To illustrate this, recently the finite temperature phase transition in the model 
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Figure 9: Pair susceptibility function of L. 



with A = 1 was studied ||11|| . In this case, the Hamilton operator can be re- written 

as 

1\ / 1\ / 1\ / 1 



^-^ r^ T "2; 

+ 4 J x -Jy-Jljl\ 

~ 4 £ ( n -T-^)(^i-^)- (3-16) 

The spin symmetry and the conservation of fermion number is obvious. Further, 
there is an on-site attraction between electrons of opposite spins like in the attractive 
Hubbard model. In figure |S] we show results for the pair susceptibility as a function 
of the spatial size L which demonstrate that the transition in this model indeed 
belongs to the K-T universality class. 

The helicity modulus is also a useful observable to establish the K-T universality 
class. It can be shown that it satisfies the finite size scaling form 

= 2 + JMT) coth (yfAtF) log(L/L (T))) , (3.17) 
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with A(T C ) = [2J j. The helicity modulus has also been measured using the im- 
proved estimator discussed earlier and the results further confirm the predictions of 
Kosterlitz and Thouless. More details can be found in Ref. UTT . 



3.5 Inclusion of a Chemical Potential 



Until now we have only considered models which can be written in terms of clus- 
ters that satisfy the three requirements necessary for the meron-cluster techniques 
to be applicable. However, in certain cases it is possible to relax these stringent 
requirements and construct models which still have no sign problem. Let us now 



demonstrate this by extending the Hubbard-type model discussed in section pT2| . We 
add a single-site term of the form 

H ' = E { U (*M - ~) (n xl ("*t + n xi - 1)} . (3.18) 

It is easy to see that U induces an on-site repulsion or attraction like in the Hubbard 
model and \i is just a chemical potential. We have already discussed that in the 
absence of this single-site term, the partition function of the model is given by 

Z= J2 2 2Nc W[b}, (3.19) 

[b] , zero— meron 

where the origin of 2 2N < can be traced to the four flips of the two identical clusters in 
the spin-up and the spin-down layer. Further, the reference configuration, obtained 
when both clusters are flipped identical to each other is guaranteed to be positive. 
It is straightforward to argue that the above single-site term modifies the factor 2 2Nc 
to 

n { 2 «p Hs 7 *] C08h Gs ± 2 exp [id 7 *] } < 3 ' 20 > 

where Sc is the size of the cluster and flc is the number of time-slices times the 
temporal winding number of the cluster. The negative sign should be taken for a 
meron-cluster. 

For the sign problem to remain solved it is important that the factor given 
in eq.( |3.20 ) remains positive. Clearly, this is the case when U < for any \i. 



Interestingly, since the size of the cluster Sc is always greater than or equal to 
fi c , it is possible to show that even when U > and \i < U/2 the factor remains 
positive. This means that we have found a repulsive Hubbard-type model in which 
it is possible to add a limited chemical potential before the sign problem kicks in. 
Although it is likely that the interesting physics of high-T c superconductivity is 
beyond this region, it appears worthwhile to study this novel model in this region of 
parameter space where the sign problem remains solved. It is likely that something 
interesting can be learned there. Finally, in the U — > oo limit one recovers the 
Heisenberg anti-ferromagnet. Anti-ferro magnetism is known to be a part of the 
high-T c phase diagram. 
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4 Efficiency of the Meron-Cluster Algorithm 



Representing the partition function in terms of the statistical mechanics of clusters 
interacting locally has the obvious advantage that numerical simulations can become 
efficient. In this section we compare the efficiency of the meron-cluster algorithm 
with that of conventional fermion algorithms. 

The efficiency of a numerical simulation method can be characterized by the 
computational cost (the computer time r) that one must invest in order to reach a 
given numerical accuracy. Obviously, the cost depends on the <i-dimensional spatial 
volume V = L d , where L is the number of lattice points per direction. In addition, 
one usually discretizes the Euclidean time-direction of extent (3 = Me into M time- 
steps of size e, such that the computational cost increases as one approaches the 
time continuum limit e — > 0. Besides these two obvious dependences, r typically 
also depends on the spatial correlation length £ via a dynamical exponent z that 
characterizes the severity of critical slowing down. Further, even if £ is small, the 
correlation length £ t /e in the time-direction becomes large in units of the temporal 
lattice spacing e when one approaches the time continuum limit. Hence, there is 
another dynamical critical exponent Zt that characterizes the corresponding slowing 
down. Altogether, the total computational cost is typically given by 

rcxL%(3/e)y[cC + c'(U^) Zt }- (4-1) 



Standard local algorithms for bosonic systems — for example, the Metropolis 
algorithm — have z, z t ~ 2. This makes it difficult to simulate systems close to 
a critical point, i.e. for large £, or close to the time continuum limit, i.e. for small 
e. Local over-relaxation algorithms can be fine-tuned to reach z,Zt ~ 1. Efficient 
non-local cluster algorithms, on the other hand, are far superior because they can 
reach z, z t ~ 0. Algorithms for bosonic systems typically require a computational 
effort that increases linearly with the space-time volume, i.e. x = d and y = 1. For 
simulations of fermionic systems this is typically not the case, even when there is 
no sign problem. Usually, one integrates out the fermions and simulates a bosonic 
theory with the fermion determinant defining a non-local effective action. The stan- 
dard fermion algorithm that is most popular in simulations of strongly correlated 
electron systems [|J computes the determinant of a large matrix of a size propor- 
tional to the spatial volume. This requires a computational effort proportional to 
the spatial volume cubed, i.e. x = 3d, even when there is no sign problem, and hence 

rcxL M (/3/e)[ce 2 + C '(6/e) 2 ]. (4.2) 

A modification of this algorithm proposed in ]2S[ has x = 2d. It should be noted that 
fermion algorithms that have originally been developed for QCD simulations have 
also been applied to strongly correlated electron systems without a sign problem. 



The Hybrid Monte Carlo algorithm |26] that was used in |2T] has 

rcxL M / 4 (/5/e) 5 / 4 [ce + c'6A], (4.3) 
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and the multi-boson algorithm 



pTfl used in |2IJ has 



rcxL d (P/eM + c'Z t /e]. (4.4) 

Here we have taken the optimistic point of view that one can, at least in principle, 
reach z, zt ~ 1 by fine-tuning these algorithms. Their computational cost has a 
much better large volume behavior than the standard algorithm, but they are not 
necessarily superior on the presently studied moderate size lattices. 

When there is a sign problem, the computational cost of the above algorithms 
increases even exponentially — not just as a power — in the space-time volume. As 
discussed before, 

r oc exp(20VAf), (4.5) 

such that formally x — y — oo. In practice, this means that systems with a severe 
sign problem can simply not be simulated. The meron-cluster algorithm allows us 
to deal with the sign problem extremely efficiently In particular, the computational 
cost increases only linearly with the space-time volume, such that x = d and y = 
1. Still, this does not necessarily mean that we can reach z ~ as in cluster 
simulations of bosonic models. This is due to the re-weighting step in the meron- 
cluster algorithm. In this step the clusters are locally reconnected such that multi- 
meron configurations are never generated. For example, a non-meron-cluster may 
be decomposed into two meron-clusters only if this does not increase the total meron 
number beyond two. Deciding if a newly generated cluster is a meron naively requires 
a computational effort proportional to £, such that z ~ 1 for the meron-cluster 
algorithm. However, with recent improvements in the implementation we have been 
able to reduce this effort to only grow with log(£). A nice feature of the meron- 
cluster algorithm is that — like other cluster algorithms — it can be implemented 
directly in the Euclidean time continuum pSfl . This completely eliminates all time 



discretization errors and all e-factors in the computational cost. The total computer 
time that is required to reach a given numerical accuracy with the meron-cluster 
algorithm is hence given by 

rcxL^cr + cUAn (4-6) 

and it is likely that z,z t rs 0. This is better than standard methods even when there 
is no sign problem, and exponentially better than any alternative method when 
there is a severe sign problem. 



5 Directions for the Future 



In this article we have outlined the essential steps that lead to a novel formulation 
of fermionic lattice theories. The essential idea behind the new approach is to 
find a cluster representation of the partition function. We have shown that in 
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this approach new models can be formulated for which the sign problem can be 
completely eliminated. Further, although the clusters themselves are non-local they 
interact locally and carry information about two-point correlations. This makes 
them especially attractive for numerical simulations. Given the difficulty in the 
numerical treatment of fermionic theories, this alternative approach turns out to be 
very useful. 

There are many interesting questions that remain unanswered in the case of 
strongly correlated electron systems. Among these the physics of competing ground 
states that may exist in high-T c materials as the doping concentration is increased is 
one of the most challenging. The model we have constructed appears to have many 
interesting features that are relevant for this type of questions. This includes anti- 
ferromagnetism and s-wave superconductivity at the two extremes of the parameter 
ranges. It would be interesting if this model shows d-wave superconductivity or 
striped phases in some intermediate range of parameters. Most importantly, the 
new Hubbard-type Hamilton operator is well suited for numerical work in a large 
region of parameter space. It would be interesting to study this model in this range. 
As a first step a mean field analysis may be useful. 

The superconductor-insulator transition is another interesting field of research 
that the present progress may shed some light on f29fl . The attractive Hubbard 



model with disorder has recently been used to study this type of quantum phase 



transition |30|. We believe that the Hamilton operator of eq. (|3.5|) with disorder is 
an alternative model for this purpose. This is possible because the sign problem 
remains solved in the presence of disorder. Further, since the algorithm for the 
new model is expected to be much more efficient, it may be possible to go to much 
smaller temperatures. 

It would be interesting to find extensions of the present work to fields like nuclear 
physics and to nano-science applications like quantum dots. It is easy to construct 
models with four layers representing spin and isospin which naturally describe the 
physics of nucleons. It would be useful to construct these explicitly so that we can 
be sure that the relevant symmetries can be maintained. Extensions of the cluster 
techniques to models in the continuum is another useful direction. Preliminary work 
shows that the non-relativistic free fermion model of eq. (|2.31| ) can be extended to 



the continuum. It would be interesting if this approach can be extended to other 
interacting non-relativistic free fermions. Such extensions can be useful in studying 
physics of many electrons enclosed within boundaries as in quantum dots. 
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